Local cluster aggregation models of explosive percolation 
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We introduce perhaps the simplest models of graph evolution with choice that demonstrate dis- 
continuous percolation transitions and can be analyzed via mathematical evolution equations. These 
models are local, in the sense that at each step of the process one edge is selected from a small set of 
potential edges sharing common vertices and added to the graph. We show that the evolution can 
be accurately described by a system of differential equations and that such models exhibit the dis- 
continuous emergence of the giant component. Yet, they also obey scaling behaviors characteristic 
of continuous transitions, with scaling exponents that differ from the classic Erdos-Renyi model. 
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The percolation phase transition on both lattices and 
networks is a subject of intense study, as it provides a 
model for the onset of large-scale connectivity in ran- 
dom media, such as resistor networks, porous rocks, for- 
est fires, and even social networks [HE]. It was recently 
shown via numerical simulation of graph evolution obey- 
ing an Achlioptas process that percolation transitions can 
be discontinuous [3 j. Starting from a graph of isolated 
nodes, at each step of the evolution, two potential edges 
are chosen uniformly at random, and using some pre-set 
criteria one edge is added to the graph and the other dis- 
carded. If the edge which minimizes the product or sum 
of the size of the two components that would be merged 
is chosen, then one can show the percolation transition 
is discontinuous. Specifically, the size of the largest com- 
ponent goes from sub-linear in system size n to a large 
fraction (bounded away from 0) of the entire network as a 
sub-linear number (n\ with I < 1) of edges are added to 
the graph. Although several recent papers have explored 
the intuition for the mechanisms behind this behavior, 
such as identifying that the evolution in the subcritical 
regime must keep larger components of similar size [4]- 
8 , there are not yet mathematical evolution equations 
describing the process. 

In contrast, more restricted Achlioptas processes evolv- 
ing under "bounded size rules" can be described mathe- 
matically. Such rules are constrained so that all compo- 
nents of size greater than some cutoff are treated equiva- 
lently. In [9 it was rigorously shown that graph evolution 
under bounded size rules can be accurately described in 
terms of differential equations. It is not known, how- 
ever, whether the restriction to bounded size rules leads 
to continuous or discontinuous percolation transitions. 

Here we introduce and analyze graph evolution models 
with choice that are both more physically motivated and 
simpler mathematically. In contrast to those in [3 , our 
models are local in the sense that the choice is constrained 
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to involve edges that share one vertex in common. Thus, 
the candidate edges span up to three components (rather 
than four as in [3]). We develop a system of differential 
equations describing the evolution of the components un- 
der the associated bounded size rules, which for the sim- 
plest model show that the system must reach a critical 
point. We implement our equations numerically and find 
they accurately predict the location of the percolation 
transition for the case of unbounded rules, and demon- 
strate via simulation that the transition for unbounded 
rules is discontinuous. 

We explicitly analyze two distinct local processes, al- 
though our approach can naturally be used to examine 
other similar processes. We call the simplest the adja- 
cent edge (AE) rule: at each step, a first vertex is cho- 
sen uniformly at random, and it must connect to one of 
two distinct additional vertices also chosen uniformly at 
random (thus both candidate edges are adjacent). Intu- 
itively, the first vertex is forced to connect to one of two 
random choices. Here we choose the edge that connects 
it to the additional vertex in the smaller component, ex- 
cept possibly in the asymptotically negligible case where 
the first vertex is in the same component as one of the 
other two (discussed in detail below). Typical evolution 
of the largest component of the graph, denoted C±, is 
shown in Fig. [I] In the bounded size rule version, the 
same rule above is applied unless both components for 
the two additional vertices have size larger than some 
bound K. In that case, we simply connect to the first of 
the two additional vertices. 

We follow the approach of Spencer and Wormald [9]. 
We start with an empty graph G of n vertices. Let xi{G) 
be the fraction of vertices in components of size i: 

x i {G) = \\{v.c{v) = i}\, (1) 

where c(v) is the size of the component containing v. 
Note that Xi(G) = m^(G), where n^(G) is the component 
density (number of components of size i divided by the 
system size n) often used in cluster aggregation models in 
the physics literature [6J [10] . For bounded size variations 
we will be interested in Xi(G) for i < K. 
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FIG. 1: Typical evolution of Ci/n for Erdos-Renyi (ER), Ad- 
jacent Edge (AE), and Triangle Rule (TR), for n = 10 6 . (Top 
inset: Example of three candidate edges for TR, and two can- 
didate edges for AE.) (Bottom Inset: A n (l/2, 0.2)/n vs n for 
AE and A n (1/2, 0.4) /n vs n for TR. Each data point is the 
average over 50 iid realizations, with error bars smaller than 
symbols. Dashed lines are A/n = 1.95n -0-323 for AE and 
A/n = 1.84n-°- 367 for TR.) 



We provide a mean-field analysis over all graph evolu- 
tions. Hence Xi becomes a function of time, Xi(t), where 
we scale so that a unit of time 1 corresponds to n edges; 
we suppress the dependence on t where the meaning is 
clear. We also use Sj = 1 — ^2 k< j %k\ that is, Sj is the 
weight of the tail of the distribution starting from j. The 
probability the first vertex is in a component of size i is 
Xi. The probability that the smaller of the two additional 



components has size j is sj 

for 1 < i < K, we have the differential equations: 



s j+i ~ 2x j s j 



xj. Hence, 



~dt 



j-\-k=i 



(2) 



This family of equations captures the distribution up to 
the bound K\ the total fraction of vertices in components 
of size larger than K is captured by sk+i = 1 — Ylf=i x %- 

To find the point where the phase transition occurs 
we consider the evolution of the second moment of the 
component sizes, which we denote by W. (Again, the 
dependence on t is implicit.) We define W = ^ J2 V c(v); 
this is the expected size of the component to which 
an arbitrary vertex belongs. We may also write W = 
YjiLi i<2n i = YnLi^i = ^2ili s i' ft nel P s notation- 



=i ' 

ally to let W* = W — Y^f=i^ x i\ nere ^* corresponds 
to the contributions to W from vertices in components 
larger than the bound K. Finally, when components of 
size j and k are merged, the change in W is equal to 



(j + k) 2 - j 2 - k 2 = 2jk. The evolution of W is: 

K K K 



dW 
~dt 



= EE 2 A'(4-4 + i) + E 2 ^^ 

j=l k=l 3=1 
K 

+ J2UW*(st-s 2 k+1 ) + 2(W*) 2 s K+1 . (3) 
k=l 



The four terms can be explained as follows. 1) Both 
components have size < K: the change in W is 2jk mul- 
tiplied by the respective probabilities that the first vertex 
has component size j and the smaller of other two com- 
ponents has size k. 2) The first vertex has component 
size j < K; the second is larger than K. (Note that the 



k=K+l 



kxk simplifies to W*, which we have used 
to simplify the expression.) 3) The first vertex has com- 
ponent size greater than K and the second does not. 4) 
All three components have size greater than K. 

As 2^=i 3 x j = W — W*, we can simplify to obtain 



dW 
~dt 



K 



= 2Wj2k(4-4+i) + 2WW*s K+1 , (4) 



k=l 



which can be used to show that this bounded size rule 
must eventually reach a critical point where W grows to 
infinity. Specifically, since > Sfc+i, 



dW 
~dt 



> 2WW*s K+1 . 



(5) 



Consider the first point where sk+i > e for some con- 
stant e > 0. (Since we keep adding edges, and K is a 
constant, it is straightforward to show that sk+i must 
eventually grow larger than a suitably small constant e.) 
At this point W* > eW (since the at least e fraction of 
W from large components must contribute at least eW of 
W's value), implying ^ > 2e 2 W 2 , from which it follows 
that W goes to infinity at some finite time. 

It is tempting (but somewhat unrigorous) to consider 
the limiting version of these eqns without the bound K: 



dW 
~dt 



= 2WJ2 S 1 > 2W - 

k=l 



(6) 



It is not immediately clear how to use Eq. [6] to similarly 
demonstrate a critical point for the unbounded case. 

Further details need to be dealt with to formalize the 
accuracy of the differential equations; here we refer the 
reader to [9 , which provides a full treatment for the case 
where two independent edges are chosen for each step. In 
particular, a key issue is that the differential equations 
fail to take into account redundant steps, where an edge 
joins two vertices that are already in the same compo- 
nent. The behavior for the x^'s is relatively straightfor- 
ward under bounded-size rules; the probability that two 
vertices chosen at random fall in the same component of 
size at most K is 0(K 2 /n), and the asymptotic effect 
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of such deviations does not affect convergence to the dif- 
ferential equations. The argument is more challenging 
for bounding the effect on W, as W's growth involves 
components of size larger than K\ however, by showing 
that the fraction of vertices in components of size k with 
high probability eventually falls geometrically with k (as 
detailed in [9]), similar bounds can be shown to hold. 

One additional benefit of considering local schemes is 
that various generalizations are entirely transparent. For 
example, the extension to d choices of neighbors of the 
first vertex instead of two for a given integer d yields 

d ^ = -ix i -i{s d i -s d i+1 ) + i Y *A4-4+i)\ ( 7 ) 

j-\-k=i 

^ = 2WJ2 k(s d k - s d k+1 ) + 2 W*4- + V (8) 

k = l 

Again, the limiting variation as K goes to infinity has 
the simpler form ^ = 2WY^k=i sf. 

The second process we study (suggested in [4j [11]) is 
the triangle rule (TR): at each step choose three distinct 
vertices uniformly at random, examine the triangle of 
three possible edges connecting the pairs of vertices, and 
select the edge that connects the two smallest compo- 
nents. Typical evolution of C\ for this process is shown 
m Fig. [1] The bounded size rule variant is that if all 
components have size above the bound K, we choose a 
random edge from the three; if two components have size 
above the bound K, we choose a random edge from the 
two adjacent to the smallest component. Using the same 
notation and analysis approach as for the AE rule, we can 
find differential equations for the bounded size variant; 

= — 2ix\ — 6ixfsi+i — 3ix 2 (l — Si) — 3iXiS 2 +1 
- 6ixiSi+i(l - + 6z Y, XjXkSh+i 

j-\-k=i;j<k 

+ix* /2 + 3z (9) 

j-\-k=i;j <k 

We briefly explain each term of Eq. [9] 1) All compo- 
nents have size z; 2z vertices lost. 2) Two components 
have size i, one has size greater than z. 3) Two com- 
ponents have size i, one has size less than z. 4) One 
component has size z, and two have size larger than z. 5) 
One component has size i, one has size less than i, one 
has size greater than z. 6) All three components have 
different sizes, and the smallest two sum to z. 7) All 
three components have size z/2. (This term only appears 
when z is even.) 8) The two largest components have 
equal size, and the smallest two sum to z. 9) The two 
smallest components have equal size and sum to z. (This 
term only appears when z is even.) Explaining, for exam- 
ple, the second term in more detail: we lose 2z vertices 
from Xi when two components have size z and one has 
size greater than z, and the probability of this is 3x 2 Si+i 
when we take into account the orderings of the choices. 



We again analyze how W = ^ c(v) changes: 

dW K ~ X K k K 

-^ = Y Y l2 jkxjX k s k+1 + Y Y Q 3 kx 3 x l 

3=1 k=j+l j=l k=j+l 

K-l K 

K 

+ W*Y SjxjSK+i + 2(W*) 2 s K+1 . (10) 
j=i 

The triangle setting lacks the pleasant form of the AE 
rule, but is suitable for calculation and fairly succinct. 
Here too, we can similarly create equations for merging 
the two largest components instead of the two smallest. 

For both the adjacent edge (AE) and triangle rule (TR) 
models, we solve the differential equations numerically 
using Euler's method in order to calculate, roughly, the 
location of the phase transition. We discretized time with 
steps of size 10 -6 . More sophisticated approaches using 
higher precision and error bounds could yield more pre- 
cise values, but the simple approach is sufficient for our 
current purposes. For the AE model, using a value of 
K = 400 led to an explosion in W occurring between 
times 0.794 and 0.795; for K = 600, the explosion oc- 
curred slightly later, between times 0.795 and 0.796. For 
the TR model, at K = 400 the explosion occurred be- 
tween times 0.847 and 0.848, and for K = 600 it occurred 
between times 0.848 and 0.849. This closely matches the 
results from direct simulation of the graph evolution pro- 
cesses discussed next. 

We establish the explosive nature of the transition for 
both the AE and TR models via numerical simulation 
of the underlying graph processes. We follow the ap- 
proach introduced in [3], while here providing a more 
formal and detailed explanation of the procedure. Let 
A n (7, A) denote the number of edges required for C\ to 
go from size C\ < [n 1 'J to size C\ > [An\ , for a system 
of n vertices. We wish to understand the asymptotic be- 
havior, lim^oo A n (7, A). If A n (7,A) increases linearly 
with rz, then the time difference spanned by the window, 
A n (7, A)/n, approaches a limiting constant greater than 
zero (the slope of A n (j : A) versus n). If, in contrast, 
A n (7, A) ex tt/ 3 with /3 < 1 (i.e., A n (j,A) is sublinear in 
n), then A n ( r y,A)/n — >• as n — >• oo. In other words C\ 
goes from size n 1 to size An in a time difference which 
approaches zero (shown for AE and TR in Fig.^a)). 

As shown in the inset to Fig[l] for the AE model we 
find that A n (0.5, A) - n 68 for all A e [0.1,0.3]. For the 
TR model we find A n (0.5, A) - n 63 for all A G [0.1, 0.4]. 
The lower bound should decrease as we access larger rz. 
The upper bound estimates the largest value of A for 
which the scaling is sublinear, denoted A c . Formally, 
A c = sup A [linin^oo A n (7, A)/n — >> 0], which is the size 
of the discontinuous jump in C\jn when viewed within 
this scaling window. 

We can bound the critical point for each process using 
the upper and lower boundaries of A n (7, A). Namely, 
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FIG. 2: (a) Measuring the lower and upper boundaries of A n (l/2, 0.2)/n for AE and A n (l/2, 0.4)/n for TR. (b) Component 
density ni ~ i -2 ' 1 shown for AE at t c (TR and PR are similar). More interesting is the rank-size component distribution (inset 
for ER and AE at t c ), showing the preponderance of large components for AE. Fitting for 50 < j < 50,000 yields Cj ~ j~ s , 
with (5 = —0.66 for ER and 5 = —0.90 for AE (TR and PR are similar but more noisy), (c) W versus t for TR. Inset shows 
W ~ (t c — t)~ a with red-line showing the best fit, attained with a — 1.13. The same red line is depicted in the main figure. 



we measure how to, the last time for which G\ < n 1 , 
and likewise how t\, the first time for which C\ > An, 
depend on n. As shown in Fig.[2|a), we find that to and t\ 
approach essentially the same limiting value denoted t c . 
Neither of these local models is as effective in delaying the 
onset of the giant component as the original Product Rule 
(PR) studied in [3] where the critical point t c ~ 0.888. 
For AE, t c « 0.796, while for TR, t c « 0.848. Likewise 
neither model is as "explosive" since A c « 0.6 for PR, 
A c « 0.3 for AE, and A c w 0.4 for TR. Other well-known 
processes have now been shown to have discontinuous 
Achlioptas process counterparts [T2HT4] . 

The discontinuous jump in the order parameter C\ is 
characteristic of first order phase transitions. Yet, we ob- 
serve critical scaling characteristic of second order tran- 
sitions. Figure [2^b) shows that n^, the scaled number of 
components of size i, behaves as ~ T r , with r = 2.1 
for both AE and TR (matching recently reported results 
for PR [6j [H US]). Figure [5|c) shows how W diverges 
at the critical point, behaving as W ~ \t — t c \~ a . We 
also see similar behavior for the size of the second largest 
component, C2 ~ \t — t c \~^ . Our numerical estimates are 



a = a « 1.13 for AE and TR, while a = /1 « 1.17 for PR. 
Note, we recover the standard Erdos-Renyi (ER) expo- 
nents (r = 5/2 and a = /i = 1) in our simulations of ER. 
Hybrid phase transitions have been previously observed 
for spin glasses p~6j [17] , constraint satisfaction problems 
(K-SAT) [18 , models of jamming in granular materials 
(see [T9U2T] and references therein), and /c-core percola- 
tion [22]. 

In summary we have introduced local models of graph 
evolution with choice that can be described by mathe- 
matical evolution equations and which exhibit discontin- 
uous percolation transitions with critical scaling behav- 
iors. Discontinuous percolation transitions are not yet 
fully understood. Local processes appear much simpler 
to describe mathematically and thus offer the potential 
for a system with a discontinuous percolation transition 
that is easier to analyze. 
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